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Abstract. We discuss the distribution of the gravitational force created by a Poissonian distribution of 
field sources (stars, galaxies,...) in different dimensions of space d. In d = 3, when the particle number 
A'^ —f +00, it is given by a Levy law called the Holtsmark distribution. It presents an algebraic tail for 
large fluctuations due to the contribution of the nearest neighbor. In d = 2, for large but finite values 
of A^, it is given by a marginal Gaussian distribution intermediate between Gaussian and Levy laws. It 
presents a Gaussian core and an algebraic tail. In d = 1, it is exactly given by the Bernouilli distribution 
(for any particle number A'^) which becomes Gaussian for ^ 1. Therefore, the dimension d = 2 is critical 
regarding the statistics of the gravitational force. We generalize these results for inhomogeneous systems 
with arbitrary power-law density proflle and arbitrary power-law force in a d-dimensional universe. 

PACS. 05.20.-y Classical statistical mechanics 



1 Introduction 

In this paper, we study the statistics of the gravitational 
force arising from a random distribution of field sources 
(stars, galaxies,...) in arbitrary dimensions of space d. This 
systematic study has interest both in astrophysics and 
probability theory. In addition, the distribution of the 
gravitational force presents numerous analogies with other 
physical systems such as Coulombian plasmas, 2D point 
vortices, dislocation systems etc. Many results have al- 
ready been obtained by Chandrasekhar [T] for the gravita- 
tional force in c? = 3 dimensions. In view of the fundamen- 
tal nature of this problem and its potential applications in 
various areas of physics and astrophysics, it is important 
to formulate the mathematical problem at a general level 
and study how the results are affected by the dimension 
of space. 

If we consider stellar systems such as globular clus- 
ters or elliptical galaxies, the problem is clearly three- 
dimensional and the gravitational force between two stars 
scales like . The structure of self-gravitating isother- 
mal and polytropic spheres has been discussed by Emden 
[2] and Chandrasekhar [3] and the thermodynamics of stel- 
lar systems has been initiated by Antonov [3] and Lynden- 
Bell & Wood [5], and developed by several authors since 
then (see the reviews of Padmanabhan [5] and Chavanis 
[7]). On the other hand, the statistics of the gravitational 
force produced by a random distribution of stars has been 
studied by Chandrasekhar [I] by analogy with the work of 
Holtsmark [5] on the distribution of the electrostatic field 
in a plasma composed of simple ions. In a series of papers. 



Chandrasekhar & von Neumann |9ll0lllll2ll3j pursued 
this work in order to obtain, from a fully stochastic the- 
ory, an expression of the diffusion coefficient of stars in a 
cluster and understand the origin of the logarithmic diver- 
gence at large scales arising in the kinetic theory of stellar 
systems (see Kandrup [H] for a review). Numerical exper- 
iments have been conducted by Ahmad & Cohen [T5] and 
more recently by Del Popolo [H] to test the predictions of 
this theory and take into account finite size effects. The 
initial theory was developed in the ease of stars but the 
same methods can also be used in cosmology assuming 
that the field sources are galaxies rather than stars |T7] . 

The standard results of Chandrasekhar [T] are valid 
for the gravitational force in three dimensions. However, 
it is important to note that some astrophysical systems 
have symmetries that lead to an effective gravitational 
interaction of lower dimensionality. 

For example, some authors have considered the gravi- 
tational interaction between infinitly elongated cylindri- 
cal filaments. In that case, the force between two fila- 
ments scales like 1/r corresponding to the gravitational 
interaction in two dimensions. The structure of polytropic 
and isothermal cylinders has been studied by Ostriker [TH] 
and the thermodynamics of gravitating rods has been de- 
veloped by Katz & Lynden-Bell [19], Aly & Perez [20] 
and Sire & Chavanis |21| . Polytropic and isothermal cylin- 
ders may have useful applications in the study of gaseous 
filaments, spiral arms and rings. Indeed, Schneider & 
Elmegreen [55] have shown that dark clouds have elon- 
gated or filamentary shapes. On the other hand, in some 
theoretical models, the spiral arms of the Galaxy are con- 
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sidered to be self-gravitating cylinders of infinite length 
|23|24j . Finally, gaseous rings occur in a variety of astro- 
nomical contexts (Saturn's ring, rings in spiral galaxies,...) 
|25|26j and infinite cylinders provide the first term of a 
natural series expansion in which one may develop the 
theory of the equilibrium of such rings. 

On the other hand, some authors have considered the 
gravitational interaction between plane-parallel sheets. In 
that case, the force between two sheets is independent 
on the distance, corresponding to the gravitational force 
in one dimension. The isothermal and polytropic distri- 
butions of such configurations have been determined by 
Spitzer and Camm [28] and systematically studied by 
Harrison & Lake [29] and Ibahez & Sigalotti [30] . On the 
other hand, their thermodynamics has been worked out by 
Katz & Lecar [31] and Sire & Chavanis [21] in the mean 
field approximation valid for N —f +oo. Interestingly, in 
the one dimensional case, Rybicki [35] has shown that the 
statistical equilibrium state can be calculated analytically 
for any N. Isothermal sheets can have application in the 
study of galactic disks, collapsing clouds, pancakes in cos- 
mology and Laplacian disk cosmogony. Indeed, in rotating 
disk systems, such as spiral and SO galaxies, the gas, dust 
and stars tend to be distributed in a symmetrical fash- 
ion about an equatorial plane. Camm |28| showed that 
the sheet model is a useful model for stellar motion in a 
direction perpendicular to the disk of a highly flattened 
galaxy. On the other hand, star-forming clouds generally 
collapse to a flattened (sheet), and sometimes filamentary 
(cylinders), configuration before fragmenting |33| . Sheet- 
like structures may also form by interstellar shocks or 
cloud collisions. Indeed, there seems to be strong evidence 
that some regions of post-shocked clouds are left near 
quasi-hydrostatic equilibrium plates (pancakes) at scales 
of galaxy formation [34] or at scales of stellar formation in 
the Galaxy [35]. Finally, plane-symmetric distributions of 
matter occur in the Saturn ring system and in the Lapla- 
cian disk cosmogony [29] . 

Apart from these various physical applications, it is 
interesting to investigate at a more academic level how 
the laws of physics, and particularly the laws of gravity, 
depend on the dimension of space d. There is indeed a 
long tradition of works in that direction |36j starting from 
a seminal paper of Ehrenfest [37| . For example, in Ref. 
[38] we have studied how the structure of relativistic white 
dwarf stars would be modified in universes with lower or 
higher dimensions and we found that the dimensions d — 2 
and d = A which surround the dimension = 3 of our 
universe are critical in some respect: white dwarf stars 
have a maximum radius in d = 2, a maximum mass in 
c? = 3 and they become unstable for c? > 4. We have also 
found that the dimensions d = 2 and d = 10 are special for 
classical isothermal spheres [21] and that the dimensions 
d = 2 and d = 9.96404372... are special for self-gravitating 
radiation in general relativity [35] ■ This type of analysis 
can shed new light on the anthropic principle and explain 
why the dimension of our universe is particular. This is a 
further motivation, in addition to the physical examples 
mentioned above, to study gravity in d dimensions. Extra 



dimensions at the microscalc also appear in theories of 
grand unification and black holes, an idea originating from 
Kaluza-Klein theory. 

On the numerical point of view, there has been con- 
siderable interest over the years in the behaviour of one 
dimensional gravitational systems, essentially by reason of 
the simplicity of these models and their relatively cheap- 
ness for numerical study. Numerical simulations are more 
easily carried out in ID than in 3-D and many early 
numerical works have considered one-dimensional self- 
gravitating systems (OGS) to study (i) the process of vio- 
lent relaxation (ii) the coUisional evolution of the system 
and its relaxation to thermal equilibrium (iii) ergodicity 
for gravitational systems. The OGS is indeed the simplest 
model for studying A^-body gravitational interactions even 
if it is not expected to capture all the features of iD in- 
teractions. Therefore, studying gravity in one and two di- 
mensions can be of interest to interprete numerical simula- 
tions. We refer to Yawn & Miller 00] for further references 
on this important topic. 

The study of the distribution of the gravitational force 
in d dimensions is also important in statistical physics and 
probability theory |41|42|43|44|45|46j because it is an in- 
teresting example of a sum of random variables where the 
Central Limit Theorem (CLT) may or may not apply de- 
pending on the dimension of space. In particular, we show 
in this paper that the dimension d = 2 is critical for the 
statistics of the gravitational field. In d = 3, the variance 
of the gravitational force produced by one star diverges 
algebraically so that the distribution of the total force is 
a particular Levy law called the Holtsmark distributior[3. 
It presents an algebraic tail which is essentially due to 
the contribution of the nearest neighbor. In d = 1, the 
variance of the gravitational force produced by one star 
is finite so that, by application of the CLT, the distribu- 
tion of the total force is Gaussian (for finite N it is ex- 
actly given by the Bernouilli distribution). In d = 2, the 
variance of the gravitational force produced by one star 
diverges logarithmically so that the distribution of the to- 
tal force is a marginal Gaussian distribution intermediate 
between Gaussian and Levy laws. It has a Gaussian core 
as if the CLT were applicable (but the variance diverges 
logarithmically with N) and an algebraic tail produced 
by the nearest neighbor as for a Levy law. Therefore, by 
changing the dimension of space, we can pass from Gaus- 
sian (d = 1) to Levy [d = 3) laws with an interesting limit 
case {d = 2). This transition has not been reported before 
in the context of gravitational dynamics and we think that 
it deserves a particular discussion. 

Finally, the systematic study of the distribution of the 
gravitational force in various dimensions of space is inter- 
esting in view of the different analogies with other physical 
systems. In d = 3, we have already mentioned the analo- 
gies between the statistics of the gravitational force cre- 

^ It is interesting to note that Chandrasekhar (1943) [1] did 
not mention the connection between the Holtsmark distribu- 
tion and Levy laws. At that time, the work of Levy (1937) [41| 
was essentially known among mathematicians and had not dif- 
fused yet in the physical and astrophysical communities. 
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atcd by stars in a galaxy and the statistics of the electro- 
static force created by a plasma composed of simple ions. 
In d = 2, the statistics of the gravitational force is similar 
to the statistics of the velocity field created by a random 
distribution of point vortices studied by Chavanis & Sire 
|47I48I49| (see also |50I51I52I53I54I55] ). to the statistics of 
the force created by dislocations [56] or to the statistics of 
the electrostatic field in a two-dimensional plasma |57| . In 
d = 1, we are not aware of any particular analogy except 
with one dimensional plasmas. 

This paper is organized as follows. In Sec.[2l we recall 
the main lines of the calculation of the distribution of the 
gravitational force in d = 3 dimensions leading the the 
Holtsmark |8llj distribution and discuss the main proper- 
tics of this distribution. In Sec. [3l we determine the dis- 
tribution of the gravitational force in d = 2 dimensions by 
adapting the results of Chavanis & Sire [33 for point vor- 
tices to the present context. In Sec. (H we determine the 
distribution of the gravitational force in d = 1 dimension. 
Finally, in Sec. [3 we generalize the results of this paper 
to the case of an inhomogeneous distribution of particles 
with arbitrary power-law density profile (or fractal distri- 
bution) and arbitrary power-law force in a d-dimensional 
universe. In Appendix [^ we give the distribution of the 
force created by the nearest neighbor in arbitrary dimen- 
sion d. For d > 2, this expression provides a good approxi- 
mation of the tail of the true distribution. Concerning the 
notations employed in this paper, wc shall use the expres- 
sion of the gravitational force in d dimensions obtained 
from the Poisson equation written as Z\<? = SdGp where 
Sd is the surface of a unit sphere in d dimensions and p is 
the density distribution. Therefore, the force (by unit of 
mass) created at the origin O by a single star located in 
r is Gm r/r^ in d = 3, Gm r/r^ in d = 2 and Gm sgn(x) 
in d = 1. With this convention, the gravitational constant 
G depends on the dimension of space (it will be denoted 
Gd in case of ambiguity). We shall also call the particles 
giving rise to the gravitational force "stars" although they 
can be filaments, sheets or other objects. 



2 The Holtsmark distribution in = 3 

The statistics of the gravitational force in d = 3 dimen- 
sions was first studied by Chandrasekhar [l] by analogy 
with the statistics of the electrostatic force in a plasma 
studied by Holtsmark [8]. Chandrasekhar computed the 
probabihty density W^(F) that a test star experiences a 
random force per unity mass F. He assumed that there are 
no correlation between the positions of the stars (Poisson 
distribution) and that the medium is infinite and homo- 
geneoufl The case of a finite uniform density distribution 



As is well-known, an infinite and homogeneous distribution 
of stars is not a steady state of a self-gravitating system. Thus, 
Chandrasekhar made a sort of "Jeans' swindle" [58]. However, 
an infinite and homogeneous distribution of masses is justi- 
fied in cosmology because the expansion of the universe has 
an effect similar (in the comoving frame) to that of a neutral- 
izing background in plasma physics [17]. On the other hand, 



has been considered by Ahmad & Cohen [T^: in that case, 
the distribution of force (F) depends on the total num- 
ber of stars N . On the other hand, the case of a system of 
non-uniform density has been treated by Kandrup |14l59j . 
His study demonstrates that the basic results are rather 
independent upon the density profile. Chandrasekhar & 
von Neumann have used their stochastic model to deter- 
mine the speed of fiuctuations T{F) [1^, the diffusion co- 
efficient of stars [To] (a calculation completed by Kandrup 
|14j ) and the spatial [ll|12j and temporal [13] correlations 
of the gravitational field. 

Let us consider a collection of N stars with mass m 
randomly distributed in a sphere of radius R with a uni- 
form density n = "iN / {AirR^) in average. The exact grav- 
itational force by unit of mass created at the center O of 
the domain is 



JV 



Gm 



Vf f. 



(1) 



In each realization, we choose at random the position of 
the stars with a uniform distribution. Since the positions 
of the individual stars fluctuate from one realization to 
the other, the value of the total force fluctuates too and 
we are interested by its distribution W(F). The problem 
then consists in determining the distribution of a sum of 
random variables. The distribution of the force created by 
one star is obtained by writing W{f)di = T(r)dr where 
t{y) = 3/(47ri?'^) denotes the density probability of find- 
ing the star in r and f = Gmr/r^ according to Eq. ([IJ. The 
Jacobian of the transformation r ^ f is readily evaluated 
leading to df = 2{Gm)-^/'^ f'^dv. Therefore, the distri- 
bution of the individual forces is given, for / > Gm/R^, 
by the pure power-law 



W{i) 



2 ^ ' 4ttR^ 



r 



-9/2 



The variance of the force created by one star 



Gm 

2 



+ 00 ^ 



(2) 



(3) 



diverges algebraically due to the behaviour at small dis- 
tances r — > (corresponding to large forces / +oo). 
Therefore, the CLT is not applicable. As we shall see, 
the distribution of the total force is a Levy law known as 
the Holtsmark distribution since it was first determined by 
Holtsmark in the context of the electric field created by a 
gas of simple ions [5] . We briefly summarize the procedure 
developed by Chandrasekhar [T] to compute the distribu- 
tion of the force. This summary is useful to compare with 
the results in other dimensions. 



in Chandrasekhar's study, since the distribution of the gravi- 
tational force is dominated by the contribution of the nearest 
neighbor, it is permissible to extend the size of the system 
to infinity without strong influence on the distribution of the 
force. As shown by Kandrup [59], only the local density of 
matter close to the star under consideration is important in 
determining the distribution of the force. 
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Since there are no correlation between the stars, the 
distribution of the gravitational force for any value of N 
can be expressed as 



N 



•' i=i \ i=i ) 



(4) 



where t(v^ = 3/(47ri?'^) governs the probability of occur- 
rence of the i-th star at position . Now, using a method 
originally due to Markov, we can express the (5-function 
appearing in Eq. ^ in terms of its Fourier transform. In 
that case, Wn{F) becomes 



with 



N 



Aiv(k) = 



47ri?3 



(5) 



(6) 



|r|=0 



where we have written f = Gmr / . We shall see that the 
distribution of the force is dominated by the contribution 
of the nearest neighbor. Therefore, we can consider the 
thermodynamic limit 



N +00, R +00, n = 



coirst. (7) 



In this limit, we obtain [T]: 
1 



with 



A(k) 



W{F) 



-nC(k) 



(27r)3 



A(k) 



(8) 



C(k) 



r+oo 

/ (1 - e*'^ '') dr. (9) 

J\r\=0 



The integral can be calculated easily [T] leading to 

4 



^(k) 



a = —{2TTGmf/^n. (10) 



Therefore, the distribution of the gravitational force is 
given by the Holtsmark distribution |8I1| : 



W{F) 



1 



+00 



1.3/2 



sm{kF)kdk. (11) 



It has the asymptotic behaviours [T]: 
2 



W{F) 



1 

3^ 



15 



1 



(27rG'm)3n2 



(F^O), 



W(F) - i(GTO)3/2„F-9/2 (F ^ +00). 
Therefore, the variance of the gravitational force 
, „ /■+°° dF 



(12) 
(13) 

(14) 
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Fig. 1. The Holtsmark distribution in c? = 3. It is compared 
with the distribution of the force due to the nearest neighbor. 
This clearly shows that the tail of the distribution is produced 
by the nearest neighbor. We have set _Fb = Gmv?^^ . 



diverges algebraically because of the contribution of large 
fields ^ 1. On the other hand the average value of the 
force is fT4l: 



{F) = AF 




(15) 



The typical force exerted upon a test particle is thus of 
the magnitude GrajD"^ which might be expected to arise 
from a few particularly nearby field stars at the interstellar 
distance D ^ ■nT^I'^ . Writing this typical force as Fg = 
Gmn2/3, Eqs. ^ and ^ give the asymptotic behaviour 
of the Holtsmark distribution for F <C Fg and F Fq 
respectively. 

It is instructive to compare the Holtsmark distribu- 
tion (fTTj) with the distribution of the force created by the 
nearest neighbor given by [1]: 



W^„.n.(F)-i(Gm)3/2„F-9/2e- 
It has the asymptotic behaviour 



(16) 



(F) ^ iG3/2m3/2nF-9/2 (F ^ +00), (17) 

which is in exact agreement with the asymptotic be- 
haviour (fT3|) of the Holtsmark distribution. Therefore, the 
highest fields are produced only by the nearest neighbor. 
By contrast, in the limit of weak forces, the two distribu- 
tions disagree: whereas the nearest neighbor distribution 
vanishes exponentially, the Holtsmark distribution tends 
to a constant value This reflects the fact that in 

the case of extremely weak forces, more than one field 
star plays a significant role. The Holtsmark distribution is 
compared to the distribution of the force due to the near- 
est neighbor in Fig. [T] and we get a good agreement for 
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sufficiently large forces F ^ Fq. The typical force due to 
the nearest neighbor is Fn.n ~ Gm/D^ ^ Gmn?/^ where 
D is the average distance between stars. It is precisely of 
the same order as the average value of the force ([T5|) due 
to all the stars. More precisely, the average value of the 
force due to the nearest neighbor is [14] : 



^ Sil=i(l(fe//)il) ^-nd since the effective force de- 
creases like 1/r'*, we find a finite value. If we write 



2/3 



— I Gmn'^'^ 



QmQGmn^/^. 



(18) 



The two resuhs (HH) and (HH) differ only by a factor 1.28. 
Therefore, large forces are due essentially to the contribu- 
tion of the nearest neighbor and the effect of all other field 
stars cancels. Thus, field stars separated by distances large 
compared with the interstellar spacing contribute negligi- 
bly to the total stochastic force. This is the reason why 
Chandrasekhar & von Neumann [101 argue that the loga- 
rithmic divergence in the diffusion coefficient of stars has 
to be cut-off at the inter-spacing distance L0. 

Chandrasekhar's approach shows that only stars close 
to the star under consideration determine the fluctua- 
tions of the gravitational field (for large forces). In fact, 
Agekyan has calculated the average contribution to 
the total magnitude of the stochastic gravitational field 
due to some field star at a distance r from the star under 
consideration. He finds that it produces an effective force 
which in good approximation can be written as [62j : 



Gm 1 
72~l + r2/yl2' 



where 

A = 



4r(|) / 3 ^ 



97r V 47r 



r\ 1/2 



-1/3 ^ Q „-i/3^ (20) 



is a characteristic length of the order of the interparticle 
distance D. For weak separations, one has /e// — > Gm/r^ 
but for large separations r ^ D, the effects of individ- 
ual stars compensate each other and the resulting force is 
reduced by a factor (r/D)'^. This corroborates the heuris- 
tic argument of Chandrasekhar [1] concerning the va- 
lidity of the two-body approximation and the effective 
cancellation of the contribution from distant field stars. 
If we were naively to assume that the average value of 
the modulus of the force is additive (|F|) — X]t=i(|f«l) 
we would find the wrong result that the average force 
is infinite: (F) = N{f) = n J^°° ^iirr^dr = +oo. 
In fact, we must sum the modulus of the effective force 



This argument is controversial because it would imply that 
the logarithmic divergence in the diffusion coefficient of charges 
in a plasma must also be cut-off at the inter-spacing distance 
while the works of Cohen et al. [60] and Balescu [61] show that 
it has to be cut-off at the (larger) Debye length. Therefore, 
for a stellar system, the logarithmic divergence in the diffusion 
coefficient of stars should be cut-off at the system's size R, or 
Jeans length, which is the gravitational analogue of the Debye 
length in plasma physics. 



(F) = iV(/e//) = n 



+ 00 



Gm 



J.2 1-1-7.2/^2 



Aiir'^dr (21) 



we obtain 
(F) = 27r2 



4F(I ) / 3 ^ 



97r V 47r 



-,1/2 



Gmn^/^ = 5.349Gmn2/3, 
(22) 



which is close to the exact result On the other hand, 
the average value of the squared force is clearly additive 
since the particles are uncorrelated 



N 



N 



ij i=l i^i 



i=l 



(23) 

Therefore, we find that the variance of the gravitational 
force is infinite 

(F^) = N{f) ^nj i^] 47rr2dr = +00 (24) 

in agreement with Eq. (fT4|l . 



(19) 3 The marginal Gaussian distribution in d 



The gravitational field produced in O by an infinite rod of 
mass per unit length fi is f — 2G^v/r^. This corresponds 
to the gravitational force f = G^mvlr^ created by a mass 
m in two dimensions provided that we make the corre- 
spondance G2 = 2G/i/m. Now, the statistics of the gravi- 
tational force in d = 2 dimensions can be directly obtained 
from the work of Chavanis & Sire [17] on the statistics 
of the velocity created by a random distribution of point 
vortices in 2D hydrodynamics (it suffices to make the cor- 
respondance 7/(27r) Gm where 7 is the circulation of a 
point vortex). There are indeed remarkable analogies be- 
tween stellar systems and 2D vortices [63]. Adapting the 
procedure of Chandrasekhar & von Neumann [TUj, Cha- 
vanis & Sire [47|48|49] have used this stochastic approach 
to obtain an estimate of the diffusion coefficient of point 
vortices when their distribution is homogeneous. This is 
another manifestation of the deep formal analogy between 
stars and galaxies. Other related works on the statistics 
of the velocity created by point vortices, including direct 
numerical simulations to test the theoretical results, have 
been performed in |50l51l52l53l54[55] . 

Let us consider a collection of N particles with mass m 
randomly distributed in a disk of radius i? with a uniform 
density n ^ NI{ttR^) in average. The force by unit of 
mass created at the center O of the domain is 



N 



Gm 



(25) 



i=l 
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The problem consists in determining the distribution of a 
smTL of random variables. The distribution of the force cre- 
ated by one star is obtained by writing W{{)di ~ r(r)c?r 
where T(r) = l/ijrE?) denotes the density probability of 
finding the star in r. Using df = {Gm)~^ f^dr, we obtain, 
for / > Gm/R, the pure power-law: 



The variance of the force created by one star 

2 f+oo 



if 



1 



f Gm 
ttR^ J \~ 



2TTrdr 



(26) 



- (27) 







diverges logarithmically due to the behaviour at small and 
large distances (corresponding to weak / and large 
/ +00 forces). Therefore, strictly speaking, the CLT is 
not applicable. However, since the divergence of the vari- 
ance is weak (logarithmic) we shall see that the distribu- 
tion of the total force is intermediate between Gaussian 
and Levy laws. The core of the distribution is Gaussian 
as if the CLT were applicable (but the variance diverges 
logarithmically with N) while the tail is algebraic, and 
produced by the nearest neighbor, as for a Levy law. 

Following the method previously exposed and consid- 
ering the thermodynamic limit 

N 

N ^ +00, R^+oo, n = — — J = const. (28) 



we obtain 



with 



A{k) 



W(F) 



iC(k) 



(27r)2 



A(k)e 



-ik-F 



(29) 



C(k) 



|r|=0 



(1 - e*-^) dr, (30) 



where we have written f = Gmr/r^. Note that we cannot 
let R +00 in the last integral since it diverges loga- 
rithmically for large r. Still, the procedure is well-defined 
mathematically if we view (|29 p - pop as an equivalent of 
Wn{F) for large N, not a true limit. The integral in (|30p 
can be calculated explicitly [47] leading to 



^(k) 



^Tr{Gm)^n. (31) 



For F < FcritiN) where F^ritiN) is defined by Eq. ([571). 
we need to consider large values of k in Eq. ([3T|) and we 
can neglect the contribution of k in the logarithm, writing 
A(k) ~ e""*^ inN _ Therefore, we get a Gaussian distribu- 
tion 



W{F) = 



1 



7i(Gm)27r2 In TV 



iF<F,„tiN)), 



(32) 



as if the CLT were applicable. However, if we were to 
extend this distribution for all values of F, we see that 
the variance of this distribution 



(33) 



diverges logarithmically with N due to cooperative effects. 
On the other hand, the average value of the force is 



1 



1/2 



(34) 



For F > Fcrit{N), we need to consider small values of k 
in Eq. (|3ip and its contribution in the logarithm becomes 



„2afe^ In fc 



In that case, we find 



crucial, so that A(k) 
after some calculation [47] that 

M^(F) = n(Gm)2F-4 {F > F,„t{N)). (35) 

Therefore, the distribution of the gravitational field in d ~ 
2 has an algebraic tail as for a Levy law. The variance of 
the gravitational force 



+00 



dF 



(36) 



diverges logarithmically due to the contribution of large 
field strengths. Comparing Eqs. (j32p and ([SS]) . we find 
that the typical force where the two regimes (Gaussian 
core and algebraic tail) connect each other is 

FcrttiN) - {nG^m^TrliiNy/^ln^^{\iiN). (37) 

For N -f-oo, Fcrit{N) +oo so, strictly speaking, the 
algebraic tail is rejected to infinity and the limit distri- 
bution W{F) is Gaussian. However, for large but finite 
values of the convergence to the limit distribution is 
so slow that the algebraic tail is always visible in practice. 
The contribution of the high field tail to the average value 
of the force is 



(F) = 2iin{Gmf 



+00 



dF 



AirnG'^Tv? 



F^^^, F^ VlnA^ln(lniV) 



1/2 



(38) 



which is smaller than the contribution ([34)) due to the core 
of the distribution. 

It is instructive to compare the marginal Gaussian dis- 
tribution with the distribution of the force created by the 
nearest neighbor given by [47j : 



^2 r-4„ ^ — 



(39) 



W„.„.(F) = n(Gm)^F-^e- 

It has the asymptotic behaviour 

M^„.„.(F) ~ n{GmfF-^ (F ^ +oo), (40) 

which is in exact agreement with the asymptotic be- 
haviour ([55)1 of the marginal Gaussian distribution. There- 
fore, the highest fields are produced only by the nearest 
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neighbor as for a Levy law. However, the Gaussian distri- 
bution of the core is created by all the particles so that the 
distribution (j39p does not provide a good approximation 
of the distribution for intermediate values of the force. The 
marginal Gaussian distribution is compared to the distri- 
bution of the force due to the nearest neighbor in Fig. [2] 
and we get a good agreement only in the tail of the dis- 
tribution. The typical force due to the nearest neighbor is 
Fn.n ~ Gm/D ^ Gm'n}/^ where D ^ rT^/'^ is the average 
distance between stars. More precisely, the average value 
of the force due to the nearest neighbor is 



{F)r 



(41) 



It is less than the average value of the force ((34|) due to all 
the stars because of the InA^ factor arising from cooper- 
ative effects. However, apart from this logarithmic term, 
they are of the same order of magnitude. This means that 
the force created by the nearest neighbor is of the same 
order as the force due to all the other particles (up to a 
logarithmic correction). This is another manifestation of 
the fact that we lie at the frontier between Gaussian and 
Levy laws. 

The present approach shows that only stars close to 
the star under consideration determine the fluctuations of 
the gravitational field (for large forces). In fact, by adapt- 
ing the calculations of Agekyan in d = 2, it is possible to 
show that the "effective" force created by a star at dis- 
tance r from the star under consideration is given in good 
approximation by [49j : 



Gm 



1 



r l + r/A' 



where 



A = (16nlniV)"^/2 



(42) 



(43) 



is of the order of the interparticle distance D. For weak 
separations, one has /e// — > Gm/r but for large separa- 
tions r ^ D, the effects of individual stars compensate 
each other and the resulting force is reduced by a factor 
(r/D). The average value of the modulus of the force can 
be obtained by summing the modulus of the effective force 
writing 

(F) = N{f,ff) = n / ———^2nrdr. (44) 



r 1 + r/A 



d=2 
ln(N)=10 



Marginal Gaussian 
distribution 




Nearest neighbor 
approximation 



10 



15 



20 



F/F„ 



Fig. 2. The marginal Gaussian distribution in d = 2. It is 
compared with the distribution of the force due to the near- 
est neighbor. This illustrates the fact that the tail of the dis- 
tribution is produced by the nearest neighbor. We have set 



4 The Gaussian distribution \n d = 1 

The gravitational field produced in O by an infinite sheet 
of mass per unit surface fiis f ~ 2t:G^, sgn(x). This corre- 
sponds to the gravitational force / = Gim sgn(a;) created 
by a mass m in one dimension provided that we make the 
correspondance Gi — 27rG/i/?7i. As we have mentioned in 
the Introduction, the one dimensional self-gravitating sys- 
tem (OGS) has been suggested as a model for the motion 
of stars perpendicular to the plane of highly flattened disk 
galaxies and it has been extensively studied in numerical 
simulations for its computational ease [40j . 

The statistics of the gravitational force in d = 1 di- 
mension created by a Poissonian distribution of stars is 
relatively straightforward. Let us consider a collection of 
N particles with mass m randomly distributed in an in- 
terval [—L,+L\ with a uniform density n = N/(2L) in 
average. The force by unit of mass created at the center 
O of the domain is 



N 

E 



ft, fi^ Gm sgn{xi), 



(47) 



This yields 



where sgn(x) = +1 if a; > and sgn(a;) = — 1 if a; < 0. 
The total force in O can be written 



(F) = ( —nG^m^n'^lnN 
16 



1/2 



(45) 



which is comparable to the exact result . On the other 
hand, the variance of the gravitational force is infinite 



= N{f) = n 



-\-oo 



2'Krdr — -foo (46) 



in agreement with Eq. p6|) . 



F ^Gm{N+- N^), 



(48) 



where iV+ is the number of stars in the interval < x < L 
and is the number of stars in the interval —L < x <0. 
Since Nj^ + = N , we can rewrite Eq. ([48|) as 



F = Gm{2N+ - N). 



(49) 



We note that, in one dimension, the gravitational force 
takes only discrete values. On the other hand, according 
to Eq. ((49|) . the probability of the fluctuation F is equal 
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to the probability of having stars in the interval < 
X < L. Since the N particles are uniformly distributed 
in the domain of size 2L, the probability that a star is 
in the interval 0<x<Lis,p~ L/{2L) = 1/2 and the 
probability that a star is in the interval —L < a; < is 
q = 1 — p = L/{2L) = 1/2. Therefore, the probability to 
have iV_|- stars in the interval < a: < L is given by the 
Bernouilli distribution: 



Nl 



N+\{N- N+y. 



N 



(50) 



The first two moments of this distribution are (A^+) = N/2 
and {{N+ - (iV+))2) N/A. Using Eq. (gS]), the distribu- 
tion of the gravitational force in d = 1 is exactly (i.e. for 
any N) given by the Bernouilli distribution 



Wn{F) 



1 



2 ' 2Gm )■ \ 2 2Gm ) ' 



(51) 



with 



{F) = 0, (F^) = NG'^m^. 



(52) 

This last result can be obtained without computation 
since the variance of the force created by one particle is 
finite and given by 



if 



(53) 



Since the particles are uncorrelated, and since (/) = 0, 

we have {F^) = E,,(/./.> = Ef=i(/f> + = 
N{p) = NG^m\ 

For 1, the CLT applies and we get the Gaussian 

distribution 



W{F) = 



1 



V2^riVGW 



^-F^/{2NG^m'^) 



(54) 



In the limit N — > +oo, the natural scaled variable is 
F/^/N. The thermodynamic limit +oo. L — > +oo 

with N/L fixed is not valid in d = 1. The result ([54]) 
can also be obtained from the Bernouilli distribution ([ST 
which becomes Gaussian in the limit of large numbers 



W{N+) 



(- 

\Nt: 



1/2 



-{N+-N/2f 



(55) 



The comparison between the Bernouilli distribution and 
the Gaussian distribution is shown in Fig. [3l Note that in 
Eqs. and the particle number iV+ and the grav- 
itational force F are treated as continuous variables so 



that the normalization conditions are W{F)dF ~ 1 
and W{N+)dN+ = 1 while in Eqs. jSO]) and ^ the 
particle number A^_|_ and the gravitational force F/Gm 
are discrete variables so that the normalization conditions 



N 



.^Wn{F) = 1 and EL=oW^Jv(iV+) 



Therefore, the relations between the discrete and the con- 
tinuous distributions are Wn{Nj^) = W{N+)dN+ with 



dN+ ~ 1 (so that Wjv(iV+) = 
W{F)dF with dF = 2GmdN+ 
2GmW{F)). 



W{N+)) and Wn{F) 
2Gm (so that Wn{F) = 



-200 - 

o 

-4001- 
-600- 



-800, 



Gaussian distribution 



Bernouilli distribution 
(exact) 



d=l 
N=1000 



200 



400 



F/F, 



600 



800 



1000 







Fig. 3. The Gaussian distribution Wo^{F) = 2GmW{F) in 
d = 1. It is compared with the exact Bernouilli distribution 
Wn[F) valid for any A^. We have set Fa — Gm and taken 
N = 1000. 



In the preceding calculations, we have assumed that 
the distribution of stars is spatially homogeneous and we 
have focused on the force at the center of the domain. 
Since the force only depends on the number of stars in 
the left and right intervals, and not on their precise distri- 
bution, the above results remain valid for any symmetrical 
distribution of the stars. If the distribution is not symmet- 
ric with respect to the point under consideration, we just 
have to compute the probability p of finding a star in the 
right interval and use the general Bernouilli formula. To 
be specific, consider an arbitrary distribution of stars with 
numerical density n{x') in the interval [Lmin, Lmax]- We 
are interested in the distribution of the gravitational force 
at X. The probability that a star is in the right interval 



[x, 



n(x') dx' , 



(56) 



and its probability to be in the left interval [Lmin,x] is 
q{x) = 1 — p{x). For example, if the stars are uniformly 
distributed in the interval [— L,L], we have p{x) = ^(1 — 
j^) and q{x) = ^(1 -|- -j). Now, the probability to have 
iV+ stars in the right interval is given by the Bernouilli 
distribution: 



WNiN+,x) 



Nl 



p{x)''+q{x) 



N+\{N -N+y. 



N~N+ 



(57) 



The first two moments of this distribution are {N+) — Np 
and {{N+ - {N+))^) ^ Npq. Using Eq. ([49]), the exact 
distribution of the gravitational force in x is 



w,iF,x)^mJpi^:^^ 



_F_V (E. _ _F^V' 

2Gm )■ \ 2 2Gm I ' 



with 



(F) =iVGm(2p- 1), 
{{F - {F)f) = ANG'^m^pq. 



(58) 

(59) 
(60) 



P.H. Chavanis: Statistics of the gravitational force 



9 



The results ([59| -(|60 |) ean be obtained without computa- 
tion. The distribution of the gravitational force created 
created by one star is W{f) = if / ^ ±Gm, W{f) = p{x) 
\i f ^Gm and W{f) = 1 - p{x) = q{x) if / = -Cm. The 
average value and the variance of the force created by one 
star are given by 

(/) = pGm + (1 - p){-Gm) = (2p - l)Gm, (61) 

if) = p(G77i)2 + (1 - p){-Gmf ^ G^m^. (62) 

Now, (F) — J^iLiifi) ^ ^(/) 8i.nd, since the parti- 
cles are uncorrelated, (F^) = J2tji.fi fj) = E!Ii(/f> + 
E.^,(.A>(/j> = N{,p)+N{N-l){f)\ Combining the pre- 
vious results, we immediately obtain Eqs. ([55)1 and (pOl) . 
Finally, for ^ 1, the CLT applies and we get the Gaus- 
sian distribution 



W{F) 



1 



^2n{{F-{F)y) 



(63) 



The result ([65]) can also be obtained from the Bernouilli 
distribution (|58p which becomes Gaussian in the limit of 
large numbers. 



5 Inhomogeneous medium and power-law 
potential in d dimensions 

In this section, we determine the distribution of the force 
created by an inhomogeneous distribution of particles in d 
dimensions. To be specific, we consider a power-law den- 
sity profile n(r) = K/r^ for r < R and n{r) = for 
r > R. We assume < p < d in order to have a decreas- 
ing density distribution that is normalizablc as r — s- 0. 
Then, K = {d - p)N / {SdR'^''^)- The uniform profile is re- 
covered for p = and K = n. For the sake of generality, 
we consider a force of the form 



f = Gm 



j.{d+a) ' 



(64) 



with d -f a — 1 > 0. The gravitational force is recovered 
for a — Q. The distribution of the force created by one 
particle is obtained by writing W{f)di = T(r)(ir where 
r(r) = {K/N)r~P. Using the transformation 

di ^ id + a- l){Gmy^n^f^^fi^ dr, (65) 
we obtain, for / > Gm/R'^^°'^^, a pure power law 

decreasing with an exponent 
7 = «i + 



d — p 



d + a-l 



(67) 



This distribution is normalizable provided that 7 > d 
which is equivalent to our previous assumptions. For the 



gravitational force [a =^ 0) and for a homogeneous distri- 
bution of stars (p = 0), one has 7 = d^ /{d — 1) for d> \. 
We note that the variance of the force created by one star 



+00 



1 



1 

^ 

j,2{d+a—l) j,p 

diverges algebraically for 



dr (X 



+00 



dr 



2a+p>2- d, 



j,d-{-2a—l-^p 

(68) 
(69) 



due to the behaviour at small distances r — > 0. In that 
case, the distribution of the (total) force is a Levy law. 
For a uniform distribution (p = 0), the criterion (|69p gives 
a > {2—d)/2 and for the gravitational case (a = 0) it gives 
p > 2 — d. For a uniform distribution and a gravitational 
force (p = a = 0), we recover the condition d > 2. In the 
following, we assume that inequality ([69]) is fulfilled. The 
critical case where (|69p is an equality, corresponding to 
a logarithmic divergence of the variance, will be treated 
specifically in Sec. 15.71 

These criteria can also be expressed in terms of the 
index 7 of the individual distribution (|66|) . For a given 
dimension of space d, we introduce the critical exponent 
7c = 2 -I- d. The variance of W{{) diverges algebraically 
(due to its behaviour for large /) when 7 < 7c. In that 
case, the distribution Vl^(F) is a Levy law. For the critical 
case 7 = 7c, the variance ofW{i) diverges logarithmically 
and W(F) is a marginal Gaussian distribution. For 7 > 7c, 
the variance of W{{) is finite and H^(F) is a Gaussian 
distribution. In the following, we assume 

d<7<7c = 2 + d (70) 

and in Sec. 15.71 we consider the critical case 7 = 7c. 



5.1 The distribution of the force 

We wish to determine the distribution of the total force 

N 

1=1 

created by the particles. Since there are no correlation 
between the particles, the distribution of the gravitational 
force for any value of N is given by 

Wn{F) = j Y{T{v,)dv,5 ( F - ^f. ) , (72) 

i=i \ i=i ) 

where Tivi) governs the probability of occurrence of the 
i-th star at position r^. Now, using the Markov method, 
we express the (5-function appearing in Eq. (|72|) in terms 
of its Fourier transform 



<5(x) 



1 



[2'kY 

Then, Wjq^) can be written 



(73) 



W'a'(F) = / A^(k)e-^'^-^dk, (74) 
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with 



N 



T(r)eJ^-^dr 



(75) 



|r|=0 



where f is given by Eq. ([64]) . Using /|^^g '''(r) dr = 1, and 
r(r) = n{r)/N, the foregoing expression is equivalent to 

A^ik) = - ^ 1^^^ (1 - e^"^ ') ^(r) dr^ . (76) 

We now consider the Umit N —>■ +oo, R +oo with 
N/R'^~P fixed. In this hmit. the distribution of the force 
can be written 



(77) 



with 



A(k) = e^^^*'^ C(k) = / n(r) (l - e'*'^) dr. 



+ 00 



(78) 



Using the transformation (j65p . we obtain 

C(k) = (Gm)^ 

J\f\=o 



(79) 



The characteristic function C(k) converges for / +oo 
if p < d and d + a-1 > 0, and for / ^ if 2a+p> 2-d. 
In other words, it converges if d < 7 < 7c. Introducing a 
spherical system of coordinates, we get 



C k = — Gm .+^-1 

ct + a — 1 



+00 



X / ^6* (sine')''"^(l-cos(fc/cos6'))/ 1 . (80) 

JO 

We have introduced the notation Cd ~ Sd/ /J^(sin0)'*~^ d6 
where Sd = 27r'*/^/I^(c?/2) represents the surface of a unit 
sphere in d dimensions. Using the identity 



we obtain 



Cd 



27r- 



r 



(81) 



(82) 



Next, setting x — kf and using the identity 

cos(a;cos6')(sin6l)'^"2d6l 
d-1 







(83) 



we can rewrite Eq. ([50]) in the form 
C(k) =afc^, 

where 

d + a — 1 



Sd 



d + a-l 



{Gm)"KB, 



(84) 



(85) 



(86) 



B = 



+00 



dx 

rH+l 



1 - r 



(87) 



The integral converges for x +00 if > 0, i.e. if d > p 
and d + a— 1>0 or, equivalently, if 7 > d. On the other 
hand, using 



1 



1 



r{i^ + i) \2i r{v + 2) \2 



y+2 



we see that the integral converges for x ^ if If < 2, i.e. 
if 2a + p > 2 — d or, equivalently, if 7 < 7c. Making an 
integration by parts and using the identity 



_d_ 

dx 



Jn{x) 



J'nix) ~ ^Jn{x) Jn+l{x) 



we can rewrite the function B in the form 



dx 

J Ax) 



(89) 



(90) 



/o X' 

It can also be expressed in terms of Gamma functions as 



1 i^/d\ r(i-f) 
2" h' \2 r{^) 



(91) 



The distribution of the force is given by Eq. ((77|) . Introduc- 
ing a spherical system of coordinates, it can be rewritten 



W{¥) 



Cd 



+00 



dk k 



d-l 



X I dO (sin ( 
/o 



\d—2 —ak —ikF cosi 

e e 



(92) 



Using identity 



wc obtain 



{2nr {2J {f 

+ 00 



e-"'' k^Ji_^{kF)dk. 



(93) 



Note that the structure of the distribution only depends 
on the scaling exponent H which takes values in the range 



Q < H < He = 2. 



(94) 
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The asymptotic behaviour of Ty(F) for small |F| can be 
obtained by expanding the Bessel function in Taylor series 
and integrating term by term. This yields 



W(F) 



1^ iw r(f + Oa^ 



(95) 



The asymptotic behaviour of M^(F) for large |F| can be 
obtained by using a method similar to the one exposed 
in Sec. II. C. of [17]. Performing the changes of variables 
z = kF and t = — cos 6 in Eq. ((92|) , we obtain 



W(F) 



Cd 



+1 



Jo 



(96) 



In this expression, t and z are real and the domains of 
integration are on the real axis — 1 < t < 1 and z > 0. 
Under this form, we cannot expand the exponential in 
power series for F — > +oo and integrate term by term 
because the integrals would diverge. The idea is to work in 
the complex plane and deform the contours of integration 
as indicated in Sec. II. C. of [17] . It is then possible to 
perform the integration on t along the semi-circle C+ of 
radius unity in the upper-half plane Im(<) > 0, and the 
integration on z along the line such that izt ~ ~y with y 
real > 0. We then obtain 

W{¥) = , ^ Rc / (1 - t^)^ dt 



r+oo 
X / e' 



Am%ryr-A y^-Uy. (97) 



We can now expand the exponential term in Taylor series 
and perform the integration on y to obtain 



(=0 •"^ + 



We now need to evaluate the integral 



r{m + d). (98) 



/ = Re 



j-Hl+d 



AHl+d 



(99) 



with t = e and 9 going from tt to 0. After straightforwad 
calculations, we can rewrite the integral (|99)l in the form 



/ = -2- 



cos 







Hl^ h - - 

2 2 2 



2 

(sin i 



Hi\e 



d0. 



(100) 



Using the identity 



cos(a6l + 6)(sin6i)'=d6' 



r{i + k) 



2'=r(i-f + |)r(i + f + |; 



■ cos 



we find that 



/ = 



.r(^) 



l+d+Hl^ 



(101) 



(102) 



Combining the previous results, the large F expansion of 
the distribution of the force can be written 

^ ' l=Q 



ttF 



l+d+HlA 
1 J 



r{-f)r{ 

We then obtain the equivalent for F +oo: 



(103) 



irF 



(104) 



Now, combining Eqs. 



and jnU), we have 



(Gm)"K^ . (105) 



Substituting this expression in Eq. (|104[1 . using Eq. 
and the identity 

2^^-'r{z)r(z + l) =^r{2z), 



(106) 



we finally obtain 



W{F) - K 



1 



d + a-1 F'i+" 



[F +oo). (107) 



Therefore, the tail of the distribution decreases with the 
exponent 



j^d + H, 



(108) 



like for the individual distribution (|66p . Furthermore, we 
shall see in Sec. 15.21 that the asymptotic behaviour of 
the distribution of the force coincides with the expression 
()118|) derived in the nearest neighbor approximation. 

The previous results can also be used for analyzing the 
stochastic gravitational fluctuations generated by a fractal 
distribution of field sources (stars or galaxies) provided 
that we make the correspondance 



p = d — df, 



(109) 
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where d/ is the fractal dimension of the distribution in a d- 
dimensional universe [64j . We also introduce the exponent 
u = d + a — 1 characterizing the power-law decay of the 
force. For the gravitational interaction i/ = d — 1. In terms 
of these quantities, the scaling exponent H can be written 

df 



H = 



The condition (|94p required to have a Levy law is 



0<df < 2v. 



(110) 



(111) 



Such a formalism can be useful in cosmology where ob- 
servations suggest that galaxies are distributed accord- 
ing to a fractal law characterized by a fractal dimension 
l<(i/<2(inac? = 3 universe) [65166] . For this range of 
values, the scaling exponent satisfies 1/2 < H < 1. 



5.2 Nearest neighbor approximation 

Let us compare these results with those obtained by mak- 
ing the nearest neighbor approximation (sec also Ap- 
pendix 1^ . For an arbitrary inhomogcneous distribution 
of particles, the distribution of the nearest neighbor is ob- 
tained from the relation 

Tn.nXr)dr = (l - _^ T„.„.(r')dr'^ n{r)Sdr''-Hr, (112) 

leading to 



r„.„.(r) = ^„(r)e-r"(-)^^-''"^'^^ (113) 

For a power-law density profile n(r) = K /r^ , we get 

T„.„.(r) = — e -^-p . (114) 

The distribution of the force due to the nearest neighbor is 
obtained from the relation iy„.„.(F)dF = r„.„.(r)rfr with 



F = Gm- 



(115) 



Using 

dF = {d + a~l){Gm)-^iTh^F^iT^ dr, (116) 
we obtain 

H.„„(F)=i^-i^!!^^^-('^+«)e-^(^r. (117) 
d + a — 1 

For F — > -|-oo, we get the asymptotic behaviour 



W„.„.(F)~/f- 



_p-{d+H) 



(118) 



d + a-1' 

which coincides with the asymptotic behaviour (|107|) of 
the exact distribution. Therefore, the tail of the distri- 
bution is dominated by the contribution of the nearest 
neighbour. We note that the moment of order b of the 
distribution (I118D is finite for b < H and its value is 



{F )n.n. 



SdK 



b/H 



(119) 



5.3 The dimension = 3 

For the ordinary dimension d = 3, we have 



H 



where 



C(k) = ak 

3-p 



H = 



47r 



2 

K 



iGm)"KB, 



i3~p)N 



(120) 

(121) 
(122) 
(123) 



As before, we assume that < p < 3, a > —2 and 2a+p > 
— 1 (or, equivalently, 3 < 7 < 5). Therefore, < H < 2. 
On the other hand, using 



-'"1/2(2;) W — sine 



(124) 



and r(3/2) = y^/2, the function B defined by Eq. dST]) 
takes the form 



B 



+00 



dx 



(x — sin x) = 



1 



H{H + 1) 



+ 00 



dx, 
(125) 



where we have used two integrations by parts to get the 
last equality. This expression can also be obtained from 
Eq. ([90]) by using the identity 



smx 



(126) 



and performing an integration by parts. The function B 
can finally be expressed in terms of the Gamma function 
as 

B = -jj-r- (127) 

2r(iJ + 2)sin(-^) ^ ' 

Using identity (|106p and the identity 

r{x)r{i - x) 



sin(7rx) 



(128) 



we can check that Eq. (|127p is consistent with Eq. (PT|) . 
The distribution of the force is then given by 



W{¥) 



1 



27r2F 



+00 



sin(fcF)fc dfc. 



(129) 



where we recall that Q < H < 2. For a = (gravity), 
we recover the situation considered by Kandrup [14] . In 
that case, H = [i — p)j2. If in addition p = (ho- 
mogeneous system), we recover the situation considered 
by Chandrasekhar [1]. In that case, H = 3/2, K = n, 
B = 4V2^/15 and a = ^{2'KGmf/'^n. We will see that 
certain results derived by Kandrup |14| contain some mis- 
takes, so it is important to reconsider this situation in de- 
tail. In addition, we treat the case of a general power-law 
potential where a can be non-zero. 
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The general properties of the distribution (|129p have we find that for — 1 < 6 < min(l,iJ) 
been derived by Chandrasekhar [67]. Although Chan- 
drasekhar considered a uniform medium, it is important 
to note that his results remain valid for a power-law den- 
sity profile; we just need to replace the index that appears 



(6 + i)r 1 



Hj r{i-b), 



' bn > 



(138) 



m his analysis hyH. Let us rewrite the results obtained by ^he link with expression ^ in the common interval of 
Chandrasekhar [67] with the present notations. When F ^ jg j^^^de by using the identity ([Ml). Finally, if we start 
IS small, a convement series expansion of the distribution directly from Eq. (24) of Chandrasekhar [67] and use the 



is given by [67 



identity p?7)l . we find that for — 3 < 6 < — 1: 



1 



+00 



In particular, 
W{¥) 



1 



:rl4 



21 



H 



1 



F 



21 



(2Z + 1)! 

(130) 



{F") = -^a"/" 

^ ' H Fi-b 



l)cos f 



-F 



h)- (139) 



27r2F \H a^l" 



1 



Finally, using the identity F{x + \) = xF{x), one easily 
checks that Eqs. and coincide in the 

{F 0). (131) whole range Q<b< H. 



Using identity (|106p . we can check that Eq. (|130p is con- 
sistent with Eq. For F ^ -\-oo. we have the series g 4 The ditnension d = 1 
expansion [67] : 



1 +°° I 



In particular, 
1 



1 



p3+lH ■ 

(132) 



For the dimension = 1 , we have 
C(k) 



where 



W{F) 



2^2 



aF{H + 2) sin [H 



1 



H = 



l-p 



2) F3 



-H ' 



Using Eqs. p22)) and p27)) . wc have 



W(F) 



p3+H ' 



(F^+oo). 



> +00). 
(133) 

(134) 



K = 



{l~p)N 



(140) 

(141) 
(142) 
(143) 



2R^-P ■ 

As before, we assume that 0<p< l,a>0 and 2a+p > 1 



This asymptotic behaviour coincide^ with the asymptotic (or. equivalently, 1< 7 < 3). Therefore, < H < 2. On 
behaviour of the distribution of the force ([TTS]) due to the ^^^^"^ ^^"^^^ '^^"^8 

nearest neighbor. Using identity (|106|) . we can check that 
Eqs. ([T^ - PM)1 are consistent with Eqs. ([TU5)) . pM)l and 
l|107p . The moments of the force are finite iff —3 < b < H. 
FoTO<b<H, we have [67] : 



COS X, 



(144) 



(ffc) = a''/"^(b + l)F{b)F (^1 
In particular, for & = 1, we get 



b 

H 



{F) = 



-F 1- — 
H 



(135) 



(136) 



and r(l/2) = ^/tt, the function B defined by Eq. ((83 
takes the form 



B = 



dx , 1 

-f^(l -cosa;) = - 



x^ 



H 



+ 00 



-dx, 



(145) 



On the other hand, if we start directly from Eq. (39) of ^^ere we have used an integration by parts to get the last 



Chandrasekhar [S7] and use the identity 



equality. This expression can also be obtained from Eq. 
([50]) by using the identity 



smcc 



dx 



2r(a)sin(^f )' 



(0 < a < 2). 

(137) 



Jl/2{x) 



sma;. 



(146) 



* Note that the expression of the asymptotic behaviour of xhe function B can finally be expressed in terms of the 
the distribution of the force obtained by Kandrup [2] is not Gamma function as 



correct. Thus, his conclusion that the asymptotic distribution 
of the force does not exactly coincide with the distribution 
produced by the nearset neighbor must be revised. 



B = 



2r(i7 + l)sin(-^) 



(147) 
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Using identities (|106p and (|128|) . wc can check that Eq. Hence, we obtain the asymptotic expansion for F 
(|147|) is consistent with Eq. ([9l|). The distribution of the 
force is then given by 



-1 I -1 

^(^) = - Ei-^y^'jnm + 1) sin {m^) 

(156) 

where we recah that < H < 2. The general proper- In particular, 
ties of the distribution p48p can be derived by adapting 

the method developed by Chandrasekhar [SJj for d = 3. W{F) ^ —ar{H + 1) sin (h—\ ^ {F +oo) 



W{F)^- e-"''" cos{kF)dk, (148) 
TT Jo 



When F is small, a convenient series expansion of W{F) ^ ' tt ^ ' V 27 F^+^ 
is obtained by expanding cos(fcF) in Taylor series and in- (157) 



Using Eqs. (|142p and (|147p . this can be rewritten 



tegrating term by term. In this manner, we obtain 
1 /9/-i-l\ 1 

In particular. This asymptotic behaviour coincides with the asymptotic 

^ / 1 \ 1 behaviour of the distribution of the force (|118[) due to the 

W{F) —F ( — j — ^, {F 0). (150) nearest neighbor. 

^ V-"/ To evaluate the asymptotic expansion of W{F) given 

Using identity we can check that Eq. ^ is con- ^y Eq. (Il5T|) for F ^ +oo, we can also integrate along 

sistent with Eq. jM]). To obtain the series expansion of imagmary axis^hus, we set z = ly with y real > 0. 

W{F) for F ^ +00, we first set z = kF and rewrite Eq. ^^^t case, Eq. ^ becomes 



M/(F) = — Rez/ e-''!^) e-^dy. (159) 



in the form 

W{F) = Rc / e-<-) e'^dz. (151) 

■^0 Expanding the first term in the integral in Taylor series. 

We now integrate in the complex plane along the line pass- '^6 get 
ing through the origin and inclined at an angle tttt to the , 

real axis instead of along the real axis itself. Thus, we set W(F) — — — Rc jgf+i v'^) ^ r(Hl + 1) (160) 

z = e^™y where y is real > in Eq. (jl5ip and we obtain t:F ^—^ l\ F^'- 



1 n+oc 

W{F) = —Rc / e-"(-F) Noting that 

T^F Jo 

,|_e-'(i-*)fy| dy. (152) 



xexp • 



Re i^'+i = - sin 



Hh 



Expanding e '^°'\ f ) in Taylor series, we have 

+ 00 



"•y F ) in Tavlor series, we have 

_ we recover Eq. (|156p . 

^/^N _ 1 \^(_jy ^g«97f From the asymptotic behaviour (|158p . it is clear that 

ttF ^ ^ *M! " F"' the moments of the force (F'') = 2/+°° W{F)F^ dF are 

finite iff — 1 < & < iJ. We can obtain an analytical expres- 

X / exp {— e^'*-^^"^^^?/! dy. (153) ^^'^ — 1 < 5 < 0. For ease of notations, we set v = —b 

Jo with < < 1. From Eq. (|148p . we have 



We rotate again the fine of integration by an angle (1 — 2 r+°° „ 

Thus, we set z = e'^^^-Ti^'iy where z is real > {F^") = - / / e""''" F'" cos{kF) dFdk. 
and we obtain ° ° 



1 a' ,r 1 

W^(^) = 'Y^^~^y ~pm We first integrate on F using the identity 

1=0 



(162) 



xeW)(i-^)fj^(^^ + l). (154) /'""S^dx^ , / , (0<a<l). 

io a;" 2r(a)cos(^^)' ^ ^ 

Now, we verify that (163) 

This yields 



Re (-z)'e'*e*(«'+i)(i-*)f = (-l)'+isin 



f+oo 

) = r^, , I e-'^^'r-ic/fc. (164) 



(155) ^ ' r{y)cos{y^) J^ 
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Expressing the integral in terms of /^-functions, we finaUy 
obtain the formula 

{F-n= f^''(^j, (0<i.<l). (165) 

HI (i^j cos(i^2 



5.5 The dimension d ~2 

For the dimension d ~ 2, we have 
C(k) = ak", 

where 

2-p 



H 



27r 



1 + a' 



l + a 



{Gm)"KB, 



K 



{2-p)N 



(166) 

(167) 

(168) 
(169) 



As before, we assume that 0<p<2,Q!>— 1 and 2a+p > 
(or, equivalently, 2 < 7 < 4). Therefore, < i? < 2. The 
function B is given by 



dx 







(170) 



Integrating by parts, the foregoing integral can be rewrit- 
ten 



1 /•+°° dx , , 



(171) 



Finally, the function B can be expressed in terms of F- 
functions under the form 



„ 1 l ^(l-f) 

2^^ffr(i + #l 



(172) 



in agreement with Eq. (|9ip . The distribution of the force 
is then given by 

1 r+00 

W{F) = — e""^" jQ{kF)kdk, (173) 

27r Jo 

where we recall that < H < 2. When F is small, a conve- 
nient series expansion of Ty(F) is obtained by expanding 
Jo{kF) is Taylor series and integrating term by term. This 
yields 



W{F) = 



1 1 

(=0 



F ( , (174) 



in agreement with Eq. (|95p . The asymptotic expansion of 
the distribution M^(F) for large F can be obtained from 
the general method developed in Sec. 15.11 leading to 



(175) 



Using Eqs. (fTBS]) . (fT72|) and the identity (fTOe]) . we obtain 
the equivalent for _F ^ +00: 



VK(F) 



(176) 



This asymptotic behaviour coincides with the asymptotic 
behaviour of the distribution of the force piSp due to the 
nearest neighbor. 

From the asymptotic behaviour (jl76p . the moments 
{F'') = J W{F)F''dF of the force exist iff -2 < 6 < iJ. We 
can obtain an analytical expression for —2 < b < —1/2. 
For convenience, we set v ~ —h with 1/2 < v < 2. Using 
Eq. (|173p and setting t = kF, we have 



+00 



Using the identity 

" Mt) 



dt 



we obtain for 1/2 < v < 2: 
1 



1 r(^^) 1 

Tlhr^ (-;^ < a < 1), (178) 

2a p f l+a \ ^2 



F-n = 



2''-^Ha''/" r(f 



(179) 



5.6 The Cauchy distribution 

We note that the characteristic function (|84p is linear 
when 

H = l, (7 = d+l). (180) 

In that case, the distribution of the force is a Cauchy law. 
This corresponds to a -I- p = 1 independently on the di- 
mension of space d. For the gravity case (a = 0), the 
Cauchy law is obtained for p ~ 1, i.e. for a fractal dimen- 
sion df = d — 1 (assuming 0? > 1). In d dimensions, the 
Cauchy law has the form 



W{F) 



\2 2 



1 



(a^ + F-^y 



(181) 



The moments (F^) exist for —d < b < 1 and they are 
given by 



(F^) = —F - + - 



1 d\F(^)Fm^) 



V2 2J F{i)r{^-±^) 
For d = 3, we obtain 

a 1 



WiF) 



with 



We have 



7r2 (a2 -fi^2)2' 



-GmK. 



W{F) 



2 + a 
1 



(F^O), 



(182) 

(183) 

(184) 
(185) 
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W(F) 



For — 3 < 6 < 1, the moments are 

^ ^ cos (^) ■ 



(F^ +oo). 



(186) 



(187) 



Using identity ^28]) . we check that Eqs. (fTST]) . (fTSS]) and 
(|182[) coincide. 

For d = 2. we obtain 



ty(F) 



with 



We have 



a 1 

2^ (a2 + ^2)3/2 ^ 

27r 



W{F) 
W(F) - 



a + 1 
1 



27ra2' 
a 



27rF3' 

For — 2 < 6 < 1, the moments are 



(F ^ 0), 
(F^ +oo). 



a" „ /I - &\ „ /6 + 2 



> V 2 
For d = 1, we obtain 

a 1 



F 



W{F) 



na^ + F^' 



with 
We have 



a = —GniK. 

a 



(188) 

(189) 

(190) 
(191) 

(192) 

(193) 
(194) 



W{F)^—, (F-^O) 
na 



W{F) 



7rF2' 

For — 1 < 6 < 1, the moments are 

(F") = 



(195) 

(F ^ +oo). (196) 



cos(&f)' 



(197) 



Using the identity ((T28)) . we check that Eqs. p97l) . p82)) 
and ()165p coincide. 



5.7 The marginal Gaussian distribution (critical case) 

We have seen that the distribution of the total force is a 
Levy law when condition (|69p . or equivalently condition 
(|7D|) . is fulfilled corresponding to < < 2. The critical 
case happens for 2a + p = 2 — d corresponding to 



H = 2, 



(7 = d + 2). 



(198) 



In that case, the variance of the force produced by a star 
diverges logarithmically. For the gravity case (a = 0) this 



corresponds to p = 2 — d. But, the condition d > 1 is re- 
quired to have a decreasing force and the condition d < 2 
must hold to have a non increasing density profile. There- 
fore, the only possibility is d = 2 and p = treated in Sec. 
[31 For other values of a, the dimension of space must lie 
in the range 1 - a < d < 2(1 - a). When H = 2, Eq. dTT]) 
remains valid with 



A(k) 



-C(k) 



C(k) / n{r) (1 - e*'^ *') dr. 

/|r|=0 

(199) 



Note that the integral defining the characteristic function 
diverges logarithmically as i? — > +oo. Therefore, Eq. ([77)1 
with Eq. (|199p must be viewed as an equivalent of the 
distribution Wn{F) for large values of R or TV, not a true 
limit. Equations ([Ml) - ([57)) are now replaced by 



C(k) = a/c^ 



where 



d + a- 1 



{GmfKB, 



(200) 



(201) 



B 



+00 



dx 



Gmk X' 



1-F 



J^_i{x) 



(202) 



Since this integral diverges logarithmically as i? -f oo, 
we can replace the term in brackets by its leading order 
expression for a; ^ using Eq. (jSH]) and we obtain 



N 



2d I G^k X 4d ^ V G2m2fc2 



(203) 



where we have used N ~ R^-p ~ R2{d+a-i) ^ 
-l-oo. Redefining 



a=- Gto) — - 

d + a - r 'Ad 

the distribution of the force can be written 



(204) 



Jo ^ 



(205) 



For F not too large (corresponding to k not too small), 
we can replace the logarithmic term by its leading contri- 
bution In A'' and we get 



W{F) = 

f-\-00 



Sd p f d 
{2tiY \2 



^~ak^ I'l^fcf J,_^(fcF)dfc. 



(206) 
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In that case, the distribution of the force is Gaussian which coincides with the distribution (jllSp due to the 

nearest neighbor for iJ = 2. We also note that the asymp- 

^^p^ _ }_ Q- i-I^N (207) totic behaviour p07p obtained for H < 2 passes to the 

{AiralnNY/^ ' hmit H ^ 2. FinaUy, the crossover between the two dis- 

tributions (|207p and (|215p occurs for a typical force 

with a variance (F^) = 2da\nN diverging like IniV. For 

F +00 (corresponding to fc ^ 0), we can replace the Fc„t{N) ^ (4a In A^)^/^ ln(ln 7V)^/^. (216) 

logarithmic term by — 21nfc and we get 



f d\ f 2\^~ It may be useful to study the physical dimensions d ~ 

^(^) ^ {2ttY^ i 2 ) \F ) 3,2,1 specifically and re-derive the previous results in a 

different manner. For the physical dimension d = 3, the 
^ / g2afe^ infc^f {kF)dk (208) expression p05p of the distribution of the force reduces to 

Jo ^ 



1 /• + °° 2, ( N \ 

Repeating the procedure of Sec. O we can rewrite the W{V) = / e ° "^Ko-^^-^k^ ) sm{kF)kdk, 

foregoing integral in the form " 



(217) 



WiF) = (^I^Re /^^(l - t')"^ dt with 



'^-(-r'"(-)e-0)'/-^y, (209) ^-fg?f- 

where we recall that the integral on t has to be performed This critical case corresponds to 2a + p = -1 or, equiva- 
in the complex plane along the semi-circle of unit radius in lently, 7 = 5. In the core of the distribution, we can make 



the upper half plane Im(t) > 0. Expanding the exponential the approximation 
term for large F ^ we obtain 



+00 ^(F) = T^/ e-'''"^^"sin(fci^)fcdfc. (219) 



1 — 2a— —r m — + 2a^^— 



^2^2 \F F'^ 



This leads to the Gaussian distribution 



W{F) = ^^^^e- ™ , (220) 



(210) (4^alniV)3/2 



Now, setting t = and integrating on C+ from = tt to ^'^^ ^ variance (F ) = 6aln7V that diverges logarithmi- 
= 0, we obtain the following results ^^Uy with N . In the tail of the distribution, we can make 



the approximation 



Re I (l-t^)'^\i'^dt = 0, (211) 1 , , 

c+ W(F) = —r- e^"'^ '"'^■sin(fcF)fcdfc. (221) 

2'K'^F , 

^ d-3 1 , Setting z = kF. this can be rewritten 

Re/ {1 - t^)— -^iUt = 0, (212) 

C+ * + T .+00 



+00 



1 



1 + ^2 / N 

M^(F) = ^— ■Im/ e^<*) ^''y^)e"zdz. (222) 



Re / ^(l-^^)^_^'^+M^ = 0, (213) 

In this expression, z is real and the domain of integration 

is on the real axis z > 0. Under this form, we cannot ex- 
pand the exponential in power series for F +00 and 
R /" fl t^)^^ '''' dt integrate term by term because the integrals would di- 

^ ' verge. However, we can modify the domain of integration 

2^-1^ 2. \ / d\ ^^'^ work in the complex plane. Then, z is viewed as a 



F ] 7^ _ 1 (214) complex variable and we can replace the domain of inte- 

r{d + 2) \ 2/ \2/ gration by the imaginary axis, i.e. z = iy with y real > 0. 

„ , . . , , . ■ -r^ JTTT7T1, i • Thus, wB havB 

Substituting these relations m hq. (|210p and using hq. 

(|204p . we obtain the equivalent for F +00: 1 r+oo _ ^ 
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Under this form, it is possible to expand the exponential 
in powers of y/F and integrate term by term. To leading 
order, we find 



W[¥) = - 



Im 



+00 



X In — 
F 



'ydy. (224) 



Since only the imaginary part of the integral matters, the 
foregoing expression reduces to 



WiF) 



+ 00 



Im / \n{i)e~yy^dy. (225) 



Then, using ln(i) = ln(e"^/2) = «7r/2 we finally obtain 

f +00 



W{F) 



2ttF^ 



which yields 



W{F) 



2ttF^ 



m 



V dy, 



3a 



ttF^' 



(226) 



(227) 



We obtain the same result if we deform the contour of 
integration in Eq. (|222p as indicated in Chandrasekhar 
[57] for the case H < 2. Substituting the value of a from 
Eq. (|218p . we obtain the equivalent 



WiF) 



K{Gm)^ 1 
2 + a F^' 



(F +00) (228) 



which is in exact agreement with the distribution of the 
force (|118p due to the nearest neighbor for iJ = 2 and 
0? = 3. We also note that the asymptotic behaviour (|134p 
obtained for H < 2 passes to the limit H ^ 2. Finally, the 
crossover between the two distributions (|220p and (|228p 
occurs for a typical force 



Fc„t{N) ^ (4alniV)i/2ln(lniV)i/2. 



(229) 



For the dimension c? = 1, the expression (|205p of the 
distribution of the force reduces to 



W{F) = 



1 



r) 



cos{kF) dk, 



with 



K(GmV 

a = 



2a 



(230) 



(231) 



This critical case corresponds to 2q: + p = 1 or, equiva- 
lently, 7 = 3. In the core of the distribution, we can make 
the approximation 



W{F) = 



1 



+00 



coa{kF)dk. (232) 



This leads to the Gaussian distribution 

1 



W{F) 



(233) 



(47ralnA^)i/2 

with a variance {F'^) = 2a\nN that diverges logarithmi- 
cally with N. In the tail of the distribution, we can make 
the approximation 

f +00 



W{F) = - 



„2afc^ In k 



cos{kF) dk. 



Setting z = kF, this can be rewritten 



W{F) = -Lrc 



+ CX3 



zHfTHf)e'^ dz 



(234) 



(235) 



To determine the asymptotic behaviour of W{F) for F — > 
+00, we deform the contour of integration and integrate 
on the imaginary axis, setting z ~ iy with y real > 0. 
Thus, we have 

W{F)^—Rer ie-2«(l')'in(^)e-2'dy. 
T^F Jo 

(236) 

We can now expand the exponential in powers of y/F and 
integrate term by term. To leading order, we find 



W{F) = 4^Re 



1 



2., J 



xlnl^ 
F 



e ^ dy. 



(237) 



Since only the real part of the integral matters, the fore- 
going expression reduces to 

2a f^°° 

WiF) ^Rc / i\n(i)e-yy'^ dy. (238) 

Then, using ln(i) = ln(e*'^/2 



it:/ 2 we finally obtain 



W{F) 



W{F) 



-a, ,2 



y dy, 



(239) 



(240) 



a r+°° 

^Jo ' 

which yields 

We obtain the same result if we deform the contour of 
integration in Eq. (|235p as indicated between Eqs. (|15ip 
and p57p for the case H < 2. Substituting the value of a 
from Eq. (|23ip . we obtain the equivalent 

K{Gm)^ 1 



W{F) 



F3' 



(241) 



which is in exact agreement with the distribution of the 
force pi8p due to the nearest neighbor for H ^ 2 and 
d = 1. We also note that the asymptotic behaviour (|158p 
obtained for if < 2 passes to the limit H ^ 2. Finally, the 
crossover between the two distributions (|233p and (I240p 
occurs for a typical force 



Fc„tiN) - (4aln7V)i/2ln(ln7V)i/2. 



(242) 



The critical case in 0? = 2 corresponding to 2a + p = or, 
equivalently 7 = 4, can be treated like in [47] . 
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6 Conclusion 

In this paper, wc have studied how the statistics of the 
gravitational force created by a random distribution of 
field sources changes with the dimension of space d. The 
dimensions d = I, d = 2 and d = 3 correspond respec- 
tively to plane-parallel (sheets), cylindrical (filaments) 
and spherical (stars) configurations. We have shown that 
the dimension d = 2 is critica|f| as it separates Gaussian 
laws (for d = 1) from Levy laws (for d > 3). This transi- 
tion may have interesting implications for the kinetic the- 
ory of self-gravitating systems (in d dimensions) since the 
distribution of the gravitational force is a key ingredient 
for the determination of the diffusion coefficient of stars 
|14|57j . Furthermore, even if our study has astrophysical 
motivations at the start, it can be of interest in probabil- 
ity theory to illustrate the differences between Gaussian 
and Levy laws. 

Note, however, that our analysis is based on several 
simplifying assumptions: 

(i) we have assumed (for the cases d = 3 and d ~ 2) 
that the number of stars +oo or equivalently that 
the system is infinite, i.e. we have considered the limit 
N,R ^ +00 with fixed n = N/V. The distribution of the 
force Wn{F) for a finite system must be computed nu- 
merically, using Eqs. ([S])-®. In d = 3 this study has been 
performed by Ahmad & Cohen [15]. It is shown that the 
convergence to the limit distribution M^(F) is quite rapid: 
for N = 2, the overall agreement with the N = oo case is 
not unreasonable, for = 50 the agreement is accurate 
to 10% and for N = 1000, the agreement is excellent. In 
d = 2, the convergence of the distributution of the gravi- 
tational force to the limit distribution (Gaussian) is very 
slow and, as we have seen, a power-law tail develops at a 
typical value of the force ([37|) which increases logarithmi- 
cally with TV. Numerical simulations exhibiting this power- 
law tail are reported (for point vortices) in |50|51|53|55j . 
In d = 1 , the exact distribution of the gravitational force 
has been obtained for any N. 

(ii) we have assumed in Sees. [2E] that the distribu- 
tion of stars is spatially homogeneous. In fact, an infinite 
and homogeneous distribution of stars is not stable and 
it clusters in dense objects (galaxies or clusters of galax- 
ies). The case of a power-law distribution of stars has been 
treated by Kandrup [M] in ci = 3 and generalized to any 
dimension and to any power-law force in Sec. [5] In prac- 
tice, since the distribution of the force is dominated by the 
nearest neighbor, the assumption of an infinite and homo- 
geneous distribution is not crucial; indeed, Kandrup [59] 
shows that the inhomogeneous case, for a smooth density 
distribution n{r), is still described by the Holtsmark dis- 
tribution (for the fluctuating force) where the density n 
is replaced by the local density n(a) at the point under 
consideration (see also Appendix A of |68] in d = 2). 

(iii) we have assumed that the positions of the stars 
are uncorrelated (Poisson distribution). This approxima- 

^ The critical nature of the dimension d = 2 has also been 
noted in [21] regarding the gravitational collapse of isothermal 
systems. 



tion may be correct in stellar dynamics where it is known 
that the two-body distribution function can be approxi- 
mated by a product of two one-body distributiont(j in a 
proper thermodynamic limit TV —f +oo [71j . However, in 
that case, the one-body distribution is spatially inhomo- 
geneous and we are led to point (ii). On the other hand, 
in cosmology, the system is statistically spatially homoge- 
neous but the particles are correlated and have the ten- 
dency to form clusters. In that case, the theoretical frame- 
work developed to determine the statistics of the gravita- 
tional force must be modified. Some interesting attempts 
to take into account spatial correlations in the position 
of the particles have been made in |72|73j . Therefore, the 
results presented in this paper can be improved in several 
directions by relaxing the above assumptions. This will be 
considered in future works. 

A last comment, suggested by the referee, may be in 
order. The limit distribution for a sum of random vari- 
ables F = X^iLi ^ classical problem in probability 
theory and there exists rigorous results and general theo- 
rems about it |41I42I43I44I45I46| . Our approach, which is 
based on the seminal work of Chandrasekhar [1] , is consis- 
tent with these general theorems. In fact, Chandrasekhar 
(1943) m and Holtsmark (1919) [8] obtained "Levy laws" 
independently from Levy (1937) [41] and other mathe- 
maticians of that time. No reference to Levy laws are made 
in the classical papers of Chandrasekhar (1943, 1948) 
[l|67j . nor in the more recent review of Kandrup (1980) 
[13] . Reciprocally, the books |42l43l44l45j do not make any 
reference to Chandrasekhar's work and derive the limit 
distributions in a more formal manner. In the present pa- 
per, we have used and extended the method introduced by 
Chandrasekhar [1] . One interest of this method is that it is 
fully explicit and amounts to the calculation of integrals. 
It is, however, restricted to pure power-law distributions 
T(f) (X /"'*' (with a cut-off at small /) while the theo- 
rems of |42|43|44|45j are more general. We have obtained 
d-dimensional generalizations of Levy laws [see Eq. ([93]) ] 
and given their main properties (usually, the problem of 
the sum of random variables is formulated in one dimen- 
sion). On the other hand, the critical case -f — jc ^ d + 2 
reported in Sees. [3| and 15.71 of the present paper (where 
the variance of the random variables (/^) diverges loga- 
rithmically) has not been treated in depth in |42|43|44|45j . 
It is usually argued that, in that case, the limit distri- 
bution M^(F) is Gaussian. This is true in a strict sense 
when N +oo. However, we have shown that for large 
but finite N, the physical distribution W^(F) has a Gaus- 
sian core and an algebraic tail. The separation between 
these two behaviours is obtained for a typical value of the 
force Fc{N) oc (IniV)^/^ ln(ln A^)^/^ which diverges with 
N. Therefore, at the limit N — > -|-oo, the power-law tail 
is rejected to infinity and only the Gaussian core remains. 
However, the convergence is so slow (logarithmic) that, in 
practice, the power-law tail is visible. This point may have 
been overlooked in [42|4 3 44 45 . 

^ Of course, for finite systems, correlations must be taken 
into account in the so-called "coUisional" regime of the dynam- 
ics, as they drive the kinetic evolution of the system |17I69I70] . 
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Finally, our study of the statistics of the gravitational 
force created by a uniform distribution of sources in d di- 
mensions illustrates the three kinds of laws that can be 
read from Fig. 1.1. of the review of Bouchaud & Georges 
[U]: (i) For d > dc = 2, the variance of the individual 
forces diverges algebraically since 7 = d^/(c? — 1) < -fc = 
2 + d [7 = 9/2 < 7c = 5 in d = 3] and the distribution of 
the total force is a d-dimcnsional Levy law ([M]) with index 
H = d/{d-l) [H = 3/2 in d = 3]. (ii) For d < 4 = 2, the 
variance of the individual forces is finite and the distribu- 
tion of the total force is Gaussian according to the CLT (in 
d = 1 it is exactly given by a Bernouilli law for all N) . (iii) 
For d = dc = 2 (critical case), the variance of the individ- 
ual forces diverges logarithmically since 7 = 7c = 4 and 
the distribution of the total force is a marginal Gaussian 
distribution. In that case, we are at the border between 
Gaussian and Levy laws (see again Fig. 1.1. of [H]). More 
generally, our results can be expressed in terms of 7 alone 
(in a space of dimension d). For d<7<7c = d-l-2, 
the variance of the individual forces diverges algebraically 
and the distribution of the total force is a d-dimensional 
Levy law ((M]) with index iJ = 7 — d (its tail decreases 
algebraically with an exponent 7 due to the nearest neigh- 
bour). For 7 > 7c, the variance of the individual forces is 
finite and the distribution of the total force is Gaussian 
according to the CLT. For 7 = 7c (critical case), the vari- 
ance of the individual forces diverges logarithmically and 
the distribution of the total force is a marginal Gaussian 
distribution. Therefore, for a fixed dimension of space d, 
the different laws correspond to d < 7 < 7c (Levy) , 7 = 7c 
(marginal) and 7 > 7c (Gaussian) . Alternatively, when we 
consider the gravitational force created by a homogeneous 
distribution of sources, 7 = d^/(d— 1) is fixed and the dif- 
ferent laws correspond to d > dc = 2 (Levy), d = dc 
(critical) and d < 2 (Gaussian). 



A The distribution of the gravitational force 
created by the nearest neighbor in d 
dimensions 

In this Appendix, we determine the distribution of the 
gravitational force in O due to the contribution of the 
nearest neighbor for a uniform distribution of stars. For 
the sake of generality, we work in a space of d dimensions. 
The probability T„.„.(r)d7' that the position of the nearest 
neighbor occurs between r and r-\-dr is equal to the prob- 
ability that no star exist interior to r times the probability 
that a star (any) exists in the shell between r and r -I- dr. 
Therefore, it satisfies an equation of the form 



This equation is readily integrated with the condition 
Tn.n.i''^) ^ Sdnr'^~^ for r — > and we find 



(245) 



From this formula, we can obtain the exact expression of 
the "average distance" D between stars: 

r + °° r+00 S^nr'^ 

D = Tn,n. {r)r dr = Sd'nr'^e 2 rfr 

Ja Jo 
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example, D = {^Y^'' 1^(4/3) in d = 3, D = 



For 

d = 2 and D = tJ- in d = 1. 

The gravitational force created in O by a star in r in 
a d-dimcnsional space is given by 



F = Gm-j. 



(247) 



We note that the variance of the force created by one star 



+ 00 
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/o ' ' ' Jo 
diverges for d > 2 due to the behaviour at small distances 
r ^ 0. In that case, the distribution of the (total) grav- 
itational force is a Levy law. For large field strengths, it 
is well-approximated by the distribution created by the 
nearest neighbor. The distribution of the force created by 
the nearest neighbor is such that Wn.n. (r)dF = r„.„. (r)dr 
where 

Tn.n.ir) = \l =ne ^ . 249 

Using 

dF = (d - l)(G7n)-'^/('*-i)F'^'/('^-i)dr, (250) 
we obtain 



^"•"■(^)='^-(d-l)F'^V(^-i) 
For F +00, we get 



W„.„.(F)~n- 



(252) 



'{d-\)FdV(d~i)- 
The distribution of the force decreases with an exponent 

d 



7 = d + 



d- 1 



(253) 



We note that the average force {F) is always finite and its 
value is 



rn.n. {r)dr ^ - ^n.n. {r')dr'^ nSdr'^'^dr, (243) „ = (^^j " p Gmn^ 



(254) 



where n is the mean density of stars. Differentiating this 
expression with respect to r, we obtain 



d_ 

dr 



Tn.n.jr) 

Sdnr"^^^ 



= ~Tn,n,{r). 



(244) 



It can be expressed in terms of the exact average distance 
between stars (12461) under the form 



(255) 
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